/**
 * Copyright 2010 - YangJiandong(chunquedong)
 * 
 * This file is part of ChunMap project
 * Licensed under the GNU LESSER GENERAL PUBLIC LICENSE(Version >=3)
 * 
 * ChunMap是开源软件。你可以自由复制、传播本项目的下载包文件，但必须保持其完整性。
 * 我们不用对使用中的风险及由此造成的损失承担任何责任。
 * 详细情况请见《ChunMap许可协议》。
 * 
 * 想了解更多有关ChunMap的信息，请访问http://code.google.com/p/chunmap/
 * 
 * 下面此段程序作者YangJianDong(chunquedong)
 * 最后修改时间2010.7.17
 */
package chunmap.model.crs.proj;

import chunmap.model.coord.Coordinate2D;
import chunmap.model.coord.Position;
import chunmap.model.coord.Transform;
import chunmap.model.crs.Projection;
import chunmap.model.crs.Spheroid;
import static java.lang.Math.*;

public class Mercator extends Projection {

	public static double maxDis = 20037508.3427892;//投影宽度的一半
    public static double r = 6378137;//地球半径

    public Mercator()
    {
        this.name = "Mercator";
    }
    
	@Override
	public Transform getReverseTransform() {
		Spheroid s=getSpheroid();
		return new MercatorRTransform(s.getA(), s.getB(), s.getE(), s.getE2(), 0, 0);
	}

	@Override
	public Transform getTransform() {
		Spheroid s=getSpheroid();
		return new MercatorTransform(s.getA(), s.getB(), s.getE(), s.getE2(), 0, 0);
	}
	
	private Spheroid getSpheroid(){
		if (getCs() != null && getCs().getSpheroid() != null) {
			return getCs().getSpheroid();
			
		} else {
			return new Spheroid("circle", r, r);
		}
	}

	private static double computeK(double a, double b, double e2, double b0) {
		double fz = pow(a, 2) / b;
		double fm = sqrt(1 + pow(e2, 2) * pow(cos(b0), 2));
		double r = (fz / fm) * cos(b0);
		return r;
	}

	class MercatorTransform implements Transform {
		private double k;
		private double l0;
		private double e;

		public MercatorTransform(double a, double b, double e, double e2,
				double b0, double l0) {
			super();
			this.k = computeK(a, b, e2, b0);
			this.l0 = l0;
			this.e = e;
		}

		@Override
		public Position convert(Position p) {
			double b = toRadians(p.getY());
			double l = toRadians(p.getX());

			double y = k * lntan(b, e);
			double x = k * (l - l0);
			
			//限制维度
            if (y > Mercator.maxDis)
            {
                y = Mercator.maxDis;
            }
            else if (y < -Mercator.maxDis)
            {
                y = -Mercator.maxDis;
            }

			return new Coordinate2D(x, y);
		}

		private double lntan(double b, double e) {
			double d = tan(PI / 4 + b / 2)
					* pow((1 - e * sin(b)) / (1 + e * sin(b)), e / 2);
			
			return log(d);
		}

	}

	class MercatorRTransform implements Transform {
		private double k;
        private double l0;
        private double e;

        public MercatorRTransform(double a, double b, double e, double e2,
                double b0, double l0)
        {
            this.k = Mercator.computeK(a, b, e2, b0);
            this.l0 = l0;
            this.e = e;
        }

        public Position convert(Position p)
        {
            double B = computeB(p.getY());
            double L = getL(p.getX());

            double b = Math.toDegrees(B);//weidu
            double l = Math.toDegrees(L);//jidu

            return new Coordinate2D(l, b);
        }

        private double computeB(double x){
            double b = getB(0, x);
            for (int i = 0; i < 20; i++)
            {
                double bb = getB(b, x);
                if (Math.abs(bb - b) < 0.01) return bb;
                b = bb;
            }
            return b;
        }

        private double getB(double b,double x)
        {
            double B=(Math.PI/2) - 2 * Math.atan(exp(b,x));
            return B;
        }

        private double exp(double b, double x)
        {
            double Exp =  Math.exp(-(x / k))
                        * Math.exp(
                                    (e / 2)* Math.log(
                                                        (1 - e * Math.sin(b))
                                                      / (1 + e * Math.sin(b))
                                                     )
                                  );
            return Exp;
        }

        private double getL(double y)
        {
            return (y / k) + l0;
        }
	}

}
